%seed
sc = parallel.pool.Constant(RandStream('Threefry'));
stream = sc.Value;
stream.Substream = 1;

Tlong = 1e5;
Xlong = zeros(Tlong, N);
shocks_indexlong = randsample(1:(T - 1), Tlong, true);

for i = 2:Tlong
    Xlong(i, :) = (Psi * Xlong(i - 1, :)')' + eps(shocks_indexlong(i), :);
end

nloop = 1e5;
results = zeros(nloop, T);

parfor iloop = 1:nloop
    stream = sc.Value;
    stream.Substream = iloop + 1;
    val = 1;

    while (val >= 1)
        init_index = randsample(1:Tlong, 1, true);
        shocks_index = randsample(1:(T - 1), T - 1, true);

        Xshort = zeros(T, N);
        Xshort(1, :) = Xlong(init_index, :);

        for i = 2:T
            Xshort(i, :) = (Psi * Xshort(i - 1, :)')' + eps(shocks_index(i - 1), :);
        end

        Psi1 = zeros(N, N);

        for i = 1:N
            regr = ols(Xshort(2:T, i), [ones(T - 1, 1), Xshort(1:T - 1, :)]);
            Psi1(i, :) = regr.beta(2:end)';
        end

        val = max(abs(eig(Psi1)));

        DR = (I_y1 + I_yspr)' * inv(eye(N) - k1x * Psi1) * X2t';
        CFX = (I_pi + I_gdp)' * Psi1 * inv(eye(N) - k1x * Psi1) * X2t';

        pdX = pxbar + CFX - DR;

        simulation_rminusgX = min(exp(-pdX));

        % reject sample if the cash flow/price ratio for GDP strip is too low (0.1%)
        if (simulation_rminusgX <= 0.001)
            val = 1;
        end

    end

    DR = (I_y1 + I_yspr)' * inv(eye(N) - k1x * Psi1) * X2t';

    CFX = (I_pi + I_gdp)' * Psi1 * inv(eye(N) - k1x * Psi1) * X2t';
    CFT = (I_pi + I_gdp + I_dt)' * Psi1 * inv(eye(N) - k1x * Psi1) * X2t';
    CFG = (I_pi + I_gdp + I_dg)' * Psi1 * inv(eye(N) - k1x * Psi1) * X2t';

    pdX = pxbar + CFX - DR;
    pdT = pxbar + CFT - DR;
    pdG = pxbar + CFG - DR;

    s1 = exp(pdT) .* taxrevgdp' - exp(pdG) .* spendgdp';

    results(iloop, :) = s1;
end

std_coeff = std(results);
